% 位移模拟器MULTI，简易演示连续体小变形的位移场
% 上一个只能设置微元体位移的版本：https://wuli.wiki/online/Strain.html
% X,Y,Z 相当于各微元顶点变形前的坐标，x,y,z相当于变形后的坐标。
% 使用连续介质力学的话术，XYZ类似Lagrange坐标和变形前的Euler坐标，而xyz相当于变形后的Euler坐标
% 你可以在弹性力学书上找到一些例子。
% CC-BY
% Gitee Repo

clc
clear

global X
global Y
global Z
global x
global y
global z

LX = 5;
LY = 1;
LZ = 1;
dX = 0.5;

% mode:'line' or 'sphere'
% 'line': 网格模式，绘制网格线
% 'sphere': 顶点/粒子模式，绘制顶点
mode = 'sphere';

% 位移场定义

% 拉伸
% ux = @(X,Y,Z) 0.1*(X+LX);
% uy = @(X,Y,Z) 0;
% uz = @(X,Y,Z) 0;

% 扭转
theta = @(X,Y,Z) pi/4*(X+LX)/(2*LX);
ux = @(X,Y,Z) 0;
uy = @(X,Y,Z) cos(theta(X,Y,Z)).*Y+sin(theta(X,Y,Z)).*Z-Y;
uz = @(X,Y,Z) -sin(theta(X,Y,Z)).*Y+cos(theta(X,Y,Z)).*Z-Z;

function h = SPHERE(x,y,z,r)
  [X,Y,Z] = sphere(20);
  h = surf(r*X+x,r*Y+y,r*Z+z);
  set(h,'EdgeColor','none')
end

function lineback(i,j,k)
    global X
    global Y
    global Z
    global x
    global y
    global z

    if i > 1
        line([X(i-1,j,k), X(i,j,k)], [Y(i-1,j,k), Y(i,j,k)], [Z(i-1,j,k), Z(i,j,k)], 'Color', 'k');
        line([x(i-1,j,k), x(i,j,k)], [y(i-1,j,k), y(i,j,k)], [z(i-1,j,k), z(i,j,k)], 'Color', 'r');
    end

    if j > 1
        line([X(i,j-1,k), X(i,j,k)], [Y(i,j-1,k), Y(i,j,k)], [Z(i,j-1,k), Z(i,j,k)], 'Color', 'k');
        line([x(i,j-1,k), x(i,j,k)], [y(i,j-1,k), y(i,j,k)], [z(i,j-1,k), z(i,j,k)], 'Color', 'r');
    end

    if k > 1
        line([X(i,j,k-1), X(i,j,k)], [Y(i,j,k-1), Y(i,j,k)], [Z(i,j,k-1), Z(i,j,k)], 'Color', 'k');
        line([x(i,j,k-1), x(i,j,k)], [y(i,j,k-1), y(i,j,k)], [z(i,j,k-1), z(i,j,k)], 'Color', 'r');
    end
end

% 生成微元体顶点网格，并计算变形后的顶点坐标
[X,Y,Z] = meshgrid(-LX:dX:LX,-LY:dX:LY,-LZ:dX:LZ);
x = X+ux(X,Y,Z);
y = Y+uy(X,Y,Z);
z = Z+uz(X,Y,Z);

% 绘制
figure();
hold on
axis equal
xlabel('x')
ylabel('y')
zlabel('z')

for i = 1:size(X,1)
    for j = 1:size(X,2)
        for k = 1:size(X,3)
            if strcmp(mode,'line')
                lineback(i,j,k);
            elseif strcmp(mode,'sphere')
                SPHERE(x(i,j,k),y(i,j,k),z(i,j,k),dX/2);
            else
                error('wrong mode')
            end

        end
    end
end

drawnow

